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I. INTRODUCTION 

Nano plasmas can now be readily produced in laser irradiated clusters, and new physical phenomena have come 
into focus experimentally as well as theoretically. Interactions between laser fields of 10 13 — 10 16 W cm -2 and clusters 
have been investigated over the last few years, see Refs. [l]-[9]. After laser interaction, extremely large absorption 
rates of nearly 100% were found in clusters, see [10], as well as X-ray radiation, see Refs. [11]— [17] . In pump-probe 
experiments by Doppncr [9] and Fennel [18] a nano plasma was generated with a first pulse in order to probe it with 
a second pulse. With the focus on optical properties, we will discuss the dynamical response function of the electrons 
in a nano plasma. For example, the absorption rate of a second laser pulse is closely connected with the dynamical 
behavior. 

Collective electronic excitations of the nano plasma are interpreted as Mie resonances of a homogeneously charged 
sphere. Absorption cross section experiments by Kresin et al. [19] show multiple resonance structures indeed. The 
effect of collective electron motion can also be seen in ultraviolet (UPS) and x-ray photoelectron spectroscopy (XPS) 
experiments, see [20], which are used to detect binding energies of core level electrons in small metal clusters [21, 22]. 
In fusion related experiments by Ditmire et al. [23] , Grillon et al. [24] , as well as Madison et al. [25] ignition processes 
are started via irradiation of deuterium clusters. Collective effects in the optical response are discussed in the context 
of metallic nanoshells by Hoflich et al. [26] as well as nanocavities by Maier et al. [27]. 

In theoretical calculations of finite systems, see Raitza et al. [28-31], a more complex resonance structure was 
found. Earlier investigations by Reinhard et al. [32] and Kull et al. [33, 34] led to comparable results. The method 
of resonance structure analysis using spherical harmonics is known from the discussion of giant dipole resonances 
of nuclei by Reinhard et al., in [35]. Quantum and semi-classical methods, see Refs. [36]-[38], respectively, were 
used to investigate the cluster excitation via laser fields. Collisional absorption processes in nano plasmas have been 
the subject of theoretical investigations in [39]. Using density functional theory (DFT) calculations, the electronic 
structures of cold clusters were analyzed by Ekardt [40], Kiimmel et al. [41], Brack et al. [42], as well as Krotscheck 
et al. [43] . The damping of collective electron oscillations was investigated by Ramunno et al. [44] emphasizing the 
importance of collisional processes beside the Landau damping. 

In this work, molecular dynamics (MD) simulations will be used to study nano plasmas in metal clusters. Clusters 
consisting of 55 up to 1000 sodium like atoms are considered after laser irradiation with intensities in the order of 
10 12 Wcm~ 2 . Properties of the nano plasma are mainly determined by the dynamics of electrons which are bound 
to the cluster but ionized from the former atoms, comparable to conduction electrons in bulk systems. As already 
shown in earlier publications, see [28], plasma parameters as known from bulk (temperature and particle density) 
but also the cluster size and net charge are justified for characterization since the electrons can be assumed to be in 
local thermal equilibrium within time scales considered here. We focus on parameter ranges where the plasma can be 
treated classically. Strong correlations are taken into account via collisions of all particles. Concepts that have been 
well established for infinite bulk systems near thermodynamic equilibrium have to be modified for applications to finite 
systems, e.g. clusters. In particular, we arc interested in the dynamical structure factor and the response function for 
such finite nano plasmas. In order to bridge from finite systems to bulk plasmas, we investigate size effects, e.g. in 
the dynamical collision frequency. First results in this direction have been reported in Refs. [28, 45, 46]. 

In Sec. II, correlation functions and their relation to optical properties of homogeneous bulk plasma are introduced 
as far as it will be of interest to extend the approach to finite systems. Expressions will also be used for comparison 
with nano plasmas in the limit of large clusters. Sec. Ill explains the restricted MD scheme for the calculation of the 
particle trajectories from which the total and bi-local current density correlation functions are determined. Symmetries 
in the correlation matrix as discussed in Sec. IV can be used for an improved statistics. The decomposition of the 
correlation matrix into eigenvectors and eigenvalues is interpreted as a decomposition into collective excitation modes. 
In Sec. V, the excitation modes will be characterized with respect to spherical harmonics. In the further analysis, we 
focus on modes with a dipole moment, which are also seen in the total current density auto-correlation function. First 
results for resonance frequencies and damping are presented. Regarding the dipole-like modes, the spatial structure 
at the selected resonance frequency will be discussed in Subsec. A and Subsec. B. Conclusion and outlook are given 
in Sec. VI. 



II. THEORY OF THE LINEAR RESPONSE OF PLASMAS IN EQUILIBRIUM 

Within linear response theory as derived by Kubo et at, see [47], the reaction of a many-particle system to weak 
external perturbations can be related to the dynamical behavior of fluctuations in thermal equilibrium. Denoting the 
equilibrium statistical operator with po, we introduce the two-time correlation function of the fluctuations SAi(t) = 
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Ai(t) — Ti{Aipo} as the Kubo scalar product 

(A i (t);A j (0))=f3 I dXTriSMWMViMPoh H,) 
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where the time dependence is given in the Hciscnbcrg picture. The indices i and j identify quantum observables. 
In particular we consider local properties so that they contain also the position r. In case of i — j, it is called an 
auto-correlation function. 

In the classical case, equilibrium two-time correlation functions can be calculated according to 

1 f T 

(A^r, ty,Aj(f*, 0)) = lim - / dr M,(r, r + t) ■ SAtf , r), (2) 

where we assumed ergodic systems - the ensemble average can be replaced by a time average. The spectrum of the 
equilibrium correlation function (A^r); Aj(f f )) u then results from Laplace transformation. 

We consider an, induced electron density fluctuation 5n e (r, t) = n e (r, t) — n e fi(r) at time t as the deviation from the 
equilibrium density distribution n e fi(f) due to an external potential U ex t{r f , t') at time t' < t, which is only dependent 
on the time difference At = t — t' and one is able to discuss its spectrum after Laplace transform. In the same way, 
the induced electric current density j e (f,t) is related to the external electric field £(? ,t'). Via Kubo's theory these 
induced quantities, 8(n e {r))^ and 6{j(r)) u , can be expressed within linear response, see [47], p. 35, as 

5(n e (r)) u = 13 j dV (5n e (f); 5h e (f<)) u U e ^,co), (3) 
S(j(f)) u = (3 ( dV <J(f); J(i*)>u, E(?, to). (4) 



The spectrum of the density fluctuation correlation (Sn e (r); 5n e (r')) UJ is related to the scalar response function. The 
current density correlation (j(f); represents in general a tensor due to the directions of the current density 

vector. The material properties of homogeneous bulk plasmas with electron density n e and inverse temperature (3 are 
only dependent on the difference of the positions Ar = r — r* . Thus, after Laplace transform of the spatial difference 
Ar, the correlations are dependent on a wave vector k. 

The optical response of a homogeneous bulk plasma will be found, i.e. discussing the dynamical structure factor, 
which is directly related to the density fluctuation correlation, as 

S(k,w) = -^{5n k ;8n k )uj- (5) 

Going from density fluctuation to current density correlation function, the dynamical structure factor can be divided 
into a static part So(k) and a dynamical part which is directly related to the longitudinal part of the current density 
correlation function 

It is of fundamental interest to describe the collective behavior of the system as response to external fields, in 
particular emission, absorption and scattering of light. In bulk systems, the wave vector and frequency dependent 
response function reads 

X (fc,o;) = -i^o^0l;il l ) w , (7) 

which can be evaluated using quantum statistical approaches such as Green function theory, see [48], or numerical 
approaches such as MD simulations, see [50]. In the long wavelength limit (k — ¥ 0), the dynamical collision frequency 
v(uj) is related to the generalized Drude formula, see [49] 

2 , ,2 



e 0 k 2 uj. 



IimY(*i,w)= t '-^-r^ , (8) 



fc->0 



(u 2 - uj 2 ^j + iuv(u) 



: 1 /2 

with the plasmon frequency co p \ = [n e e 2 / (m e e 0 )] . In the classical case the current density correlation function 
has been extensively discussed including the long wavelength limit k — > 0 applying MD simulations and perturbative 
approaches. Exemplarily, one can refer to [50]. 
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The state of a homogeneous one-component plasma in thermodynamic equilibrium is characterized by the non- 
ideality parameter T — e 2 (Aim e /?>) 1 l 3 (Air£okB,T e )~ 1 and the degeneracy parameter 0 = 2m e kBT e h~ 2 (3ir 2 n e )~ 2 / 3 . 
Considering the response function x(fc, u>) in the long wavelength limit, a sharp peak arises at the plasmon frequency 
u! p \. For finite wavelengths the resonance is shifted and can be approximated by the so called Gross-Bohm plas- 
mon dispersion for small wave numbers fc, see [51, 52], oj(k) w w p \ + 3fc 2 / k 2 + ... with the Debye screening length 
kT 1 = [n c e 2 1 (eofcsTe)] -1 / 2 . This relation has recently been revisited by [53] with respect to the relevance of collisions. 
The general behavior of the response function x(fc,w) in the long wavelength limit can be understood with the help 
of the dynamical collision frequency, using Eq.(8), see [53] and [54]. In the two-component plasma, a phonon mode 
can arise in addition to the plasmon excitations. 

The response function \(k,u>) and the related dynamical structure factor S(k,u>) as well as the optical properties 
have been intensively investigated for electron-ion bulk systems, see [53] and [55] , and will not be further reported here. 
In this work, the inhomogeneous case of finite clusters in local thermal equilibrium will be discussed. The response of 
inhomogeneous systems is not only dependent on the difference of the positions, but on r and r* separately. Therefore, 
the bi-local current density correlation function (j(r); j(r > )) UJ can not be diagonalized by spatial Fourier transform. 
Instead of plane waves, other mode functions will be found below to describe the collective excitation of electrons. 

III. MD SIMULATIONS OF FINITE PLASMAS 

Finite plasma systems have been investigated using the restricted molecular dynamics (RMD) simulations based on 
a MD scheme by Suraud et al. [56] . A two-component system of singly charged ions and electrons will be described 
using an error function pseudo potential for the interaction of particles i and j 

^«>=isHt)' <»> 

where Zj is the charge of the ith particle. The Coulomb interaction is modified at short distances, assuming a 
Gaussian electron distribution that takes quantum effects into account. Considering a sodium like system, the potential 
parameter A = 0.318 nm was chosen in order to reproduce the ionization energy of Ip = V e i(r — > 0) = —5.1 eV for 
solid sodium. 

The velocity Verlet algorithm [57] was applied to solve the classical equations of motion for electrons and ions 
numerically. This method takes into account the conservation of the total energy of the whole finite system, as long 
as there is no external potential. To follow the fast electron dynamics, time steps of 0.01 fs were taken to calculate 
the time evolution. Contrary to bulk MD simulations no periodic boundary conditions are applied. 

Icosahedral arrangements of 55, 147, and 309 ions, see [28], were considered as initial configuration for the ion 
positions. For these nearly spherically, homogeneously distributed ions, the ion density typical for solid sodium is 
given by an ionic next neighbor distance of do = 0.212 nm. In addition, randomly distributed ion configurations 
within a given sphere were considered for comparison and the number of ions was increased up to 1000 particles. 
Starting with a neutral cluster, the electrons have been positioned nearby the ions with small, randomly distributed 
deviations from the ion positions. 

To simulate experiments where clusters are excited by short pulse lasers, MD simulations are performed under 
the influence of an electric field, assuming a Gaussian shape and pulse duration of about 100 fs. Due to the largely 
increased kinetic energy of the electrons, ionization processes occur. After the laser field is switched off, the ionization 
degree of the cluster is determined by the number of electrons found outside the cluster radius with positive total 
energy, so that they can escape from the cluster. Due to ion excitation on larger time scales, a slow expansion of the 
cluster begins. 

Considering the single-time properties, it was found in [28] that already local thermodynamic equilibrium (LTE) 
is established within fs immediately after the electron heating. In particular, at each time step, the momentum 
distribution of electrons is well described by a Maxwell distribution, and the spatial density profile agrees with a 
Boltzmann distribution with respect to the average potential that is determined by the actual ion configuration and 
the self-consistent electronic mean field. The fact that electrons relax to LTE within a sub fs time interval, where the 
ion configuration remains nearly unchanged, enables us to separate the electron dynamics from the ion dynamics. 

Subsequently, the dynamical properties of the electron subsystem can be calculated for a frozen ionic configuration 
thus referring to a specific time. This is considered as an adiabatic approximation to the true dynamical properties 
of the electron subsystem which have to take into account the slow change in the ion configuration. More rigorously, 
non-stationary time dependent correlation functions have to be treated for the full charged particle system. 

Using the restricted MD (RMD) simulations scheme as introduced in [28], the ions are kept fixed acting as external 
trap potential. Starting from an initial state, the many-electron trajectory {fi(t),pi(t)} is calculated, solving the 
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classical equations of motion of the electrons. From this, all further physical properties of the electron subsystem inside 
the cluster are determined. Within RMD simulations, we consider no temporal variation of the plasma parameters 
that are determined by the frozen ion distribution, and the degree of ionization. A long-time run can be performed 
in order to replace the ensemble average by a temporal average. This has been successfully done for the single-time 
properties such as the momentum distribution and the density profile, see [28] and will now be applied to the two-time 
correlation functions. 

Using classical MD simulation techniques, the results are valid for non-degenerate plasmas. This restricts the 
temperature range to T > 1 cV where our simulations can be compared with realistic sodium clusters. Arbitrary 
values for the plasma parameter can be treated since we are not confined to the weak coupling limit as , e.g., in 
perturbation theory. 

In our RMD calculations, we start from a homogeneous ion configuration (icosahedral or randomly distributed) 
inside the cluster at fixed ion density. A thermostat was used to heat the electrons. At each time step, a heating rate 
was calculated via comparison of the electron temperature T e , calculated from the kinetic energy, and the intended 
final temperature T a i m . Hot electrons are emitted so that the cluster becomes ionized. Evaluating the trajectories 
of electrons, sufficient time of about 200 fs has to be allowed before a stationary ionization degree is established and 
thermal equilibrium is achieved, characterized by stationary distributions in position and momentum space. 

Using the trajectories of all N e electrons obtained from the RMD simulations scheme, the local current density 
j e (r,t) at position r was calculated for each time step t 

1 N " 

lm = A hm o -^X>w wnim do) 

which is the sum over all electron momenta pi inside a small volume AVp at position r, where SAV?(ri(t)) = 1, and 
^Ay ? (n(t)) = 0 for electrons found outside AVp. The size of the volume determines the spatial resolution of the local 
current density j e (r,t). However, it must be taken sufficiently large to reduce statistical fluctuations. 

The bi-local correlation tensor of the spatially resolved current density is calculated according to Eq.(2) as 

1 Nt 

(Uf,0yJe(^,t)) = ^Y,J e (r f ,i-T)®j e (r',i-T + t), (11) 

T 1=1 

with typical values of N T ~ 10 9 and r ~ 0.1 fs. Its Laplace transformation reads 

(Je(r); j e tf))» = dte iw * (j e (f,0);j e (7,t)) . (12) 

In the following, we restrict ourselves to the diagonal components (je',je)w of this tensor, where only parallel compo- 
nents of the current density vectors are correlated as already introduced in Sec. I. As it will be shown in the following 
sections, this bi-local current density correlation is important to understand the excitation modes of nano plasmas. 
The non-diagonal components of the bi-local correlation tensor are small in comparison to the diagonal components. 
Beside the bi-local current density correlation function considered here, the bi-local density fluctuation correlation 
(Sn(r),5n(r')) u as well as the bi-local force correlation {F(r) 7 F(r')) u are useful quantities in the context of optical 
properties. These correlations can be evaluated from the trajectory in a similar way and are related to the bi-local 
current density correlation, but will not discussed here any further. 

Because of the spherical symmetry of the cluster geometry during excitation and expansion, the volume is divided 
into cells AVp according to equidistant intervals of spherical coordinates, i.e. the distance to the center of the cluster 
r, the inclination angle 9 as well as the azimuthal angle <j>. The cluster radius i?i is given by the rms radius of 
ions according to Rf = 5/3(r 2 ). Because the cluster is positively charged, all bound electrons will move inside this 
cluster radius. The cluster will be divided into a finite number of sections, which are numbered by a single counter 
a = N^Nq (k — 1) + (j — 1) + i with three independent counters accordingly to the three coordinates: i = l..N^, 

j = l..Ng and k = l..N r . With respect to Eq.(ll) the bi-local correlation matrix D a<a i(t) = (j^\r a , 0); j^\r a >, t)^j 

for the spatially resolved cluster and its Laplace transform D a , a >(uj) = / 0 °° dt e lw * D a ^ a i{t) have been calculated. 

The total current density auto correlation function (ACF) can be calculated from the trajectories as well as from 
the frequency spectrum of the current density correlation matrix 

0M>« = ^V^l^Vi^A^y,^ (13) 



cl 



a,a' 



() 



using the cluster volume V c \ — ^Rf and the individual cell volumes 

rRik/Nt 



d(j) / d(9 
27v(i-l)/N 4 , J-K{j-X)/Ng 
3 

(3fc 2 - 3k + 1) ( cos 



dr~ 2 



r sin t 



2?r 



3JV. 



Ri(fe-l)/JV r 
7T 



AT* 



(i-i) 



cos 



(14) 



IV. FROM BI-LOCAL CORRELATION FUNCTION TO ELECTRON MODES 

In the following, we discuss calculations for the current density ACF Eq.(13) and the bi-local current density 
correlation spectrum D a ^ a i{oj). Exemplarily, we present results for the Na309 cluster at electron temperature T e = 1 
eV, cluster charge Z = 16 and ionic density n; = 28.0- 10 21 cm -3 . Calculations of other cluster sizes will be presented 
in the following sections. The real part of the current density ACF Re(je ;je )u i s shown in Fig. 1. Three maxima are 
obtained. To investigate the origin of the maxima as collective excitations of the nano plasma, the bi-local current 
density correlation matrix was calculated as well. 



\ 0,01 



u 0,001 




0,0001 - 



CO in fs 



FIG. 1. Frequency spectrum of the real part of the current density ACF of a Na3og cluster at electron temperature T e = 1 eV, 
cluster charge Z = 16 and ionic density rii = 2.80 • 10 22 cm -3 . 



The following spatial symmetries in the matrix D aja i(ui) were found 

Di.j,k\ i',j',k'{u) = Di* t j t f,\ j',k'(^), 

Di,j,k\ = D i,Ng-j+l,k\ V ,N e -j' +l,fc'( w ); 



(15) 
(16) 
(17) 



where i* = \ i — i'\ or i* = — \i — i'\, if \i — i'\ > N$/2. In our case, the N% elements of the full matrix can 

be reduced to Aj n d = ^ Nr Ne+ p Nr Nc> (jV^ + 7V0mod2) independent elements due to the symmetries Eq.(??), thus 
improving statistics via averaging equal elements. 

Because of the different size of the section volume there are large variations in the statistical errors of the bi-local 
current density correlation matrix elements. The statistical errors are of comparable size in the bi local current 
correlation function 



for which the following eigenproblem was solved, 



(18) 



(19) 
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Thus, the matrix is decomposed into eigenvectors $^0(0;) = 4 , A i,i I j,fe(w) as well as their eigenvalues K^(uS) at each 
frequency. The eigenvectors represent the spatial structure of the mode (W^.j j.fc(w) — > \E' M (r, w)). The orthonormality 
condition 

Jd 3 r^^(f,Lu) 1y {r,w) = . (20) 

holds. 

For two selected frequencies, the 10 strongest eigenvalues of the Na3og cluster are shown in Fig. 2 at u> = 4.42 fs -1 
(black) a resonance frequency was found, with one outstanding, leading eigenvalue. At off-resonant frequencies, i.e. 
uj = 5.50 fs -1 (red), all eigenvalues are of the same order of magnitude. 




123456789 10 



FIG. 2. Eigenvalues of the Na3og cluster sorted by size for the resonant case at cj = 4.42 fs 1 (black) and a non-resonant case 
at lo = 5.50 fs -1 (red), for parameters as in Fig. 1. 

In Fig. 3 (left) 10 eigenvalues K^(uj) of the Na3og cluster, as they result from Eq.(19), are shown in dependence 
of frequency. The are colored according to their strength and numbered ascending with descending strength. In 
the shown frequency range, modes if M (w) with well defined maxima are found. Spatial oscillation structure can be 
identified by analyzing the eigenvectors. In Fig. 3 (right), the spectra of eigenvalues are sorted in an alternative 
way, according to the same spatial oscillation structure of the eigenvector, which is obtained over the whole frequency 
range. 

It is possible to separate resonance structures of different modes (distinguished by color). The black solid mode is 
the strongest. The resonance frequencies are also found in the total current density ACF (indicated via horizontal 
blue dashed lines). Resonances in the total current density ACF, shown in Fig. 1, are only possible in the case of 
non-zero total current, which is caused by a dipole-like oscillation. Therefore, resonances which are seen in the total 
current density ACF are oscillation modes with a dipole moment. After characterization of the resonance structures, 
the dipole-like resonances will be investigated in greater detail. 



V. ANALYSIS OF THE COLLECTIVE MODES 



The decomposition of the locally resolved current correlation matrix into eigenvalues K^(lo), as shown in Fig. 
3, gives a complex set of resonance structures. To analyze the spatial oscillation structure, a spherical Fourier 
decomposition of the eigenvectors into the spherical Bessel function ji(k n jr) and spherical harmonics Y^ rn (9,(f>) was 
performed according to 

N n JV, I 

y^r,u) = J2Yl E S nJ,M N n,l3l(kn,ir)Y hl n{e,<!>), (21) 

n— 1 I— 0 m— — l 




where S n i m (u) is the spherical Fourier component with ordinal number n,l,m. The normalization factor N n ; as 
well as the wave number k n ^ are chosen in the way that the eigenfunction has a root at the cluster surface. 

In Fig. 3 (right) the four leading eigenvalue spectra were identified by the pair of the leading ordinal numbers I, m. 
The angular part of the eigenfunction to the spherical harmonics Yi im (0, (j>)- The dipole-like resonance structure is a 
mix of the same two spherical harmonics at all frequencies. As a consequence, one is able to point out similar properties 
of slightly different of eigenvectors ^(r, w) of the same resonance structure at different frequencies, represented by 
the pair of leading ordinal numbers I, m. 

For example, the dipole-like mode, represented via solid black lines in Fig. 3 (right), with resonances at the same 
frequencies as in the total current density ACF (indicated via horizontal blue dashed lines), is characterized by the 
mix of the spherical harmonic functions Y 0 fi(8, </>) and Y 2 fi(0, <j>). For the Na3og cluster, one can find three resonance 
frequencies. To draw a complete picture, the decomposition of the current matrix into eigenvectors and eigenvalues 
was done for a larger cluster consisting of 1000 ions. There, four pronounced dipole-like resonances were found. The 
spatial structures of the leading dipole-like mode of the current density j(f) ~ AV(r) * s snown f° r the Naiooo cluster 
in Fig. 4 at the resonance frequencies. The spatial mode structure of the current density j^r) ~ is shown. 



cj = 4.80 fs' 1 u = 6.18 fs" 1 uj = 9.15 fs' 1 uj = 8.23 fs" 1 




x/r d x/r d x/r d x/r cl 



FIG. 4. Selected eigenvectors of the dipol-like mode in the Naiooo clusters for same parameters as in Fig. 1 but Z — 19. 

At the resonance frequency ur = 4.80 fs, the electrons are oscillating with current densities j e {v) depending on the 
distance r, which is related to the density profile n e (r) of the electrons. As shown in Fig. 4 on the left hand side, all 
electrons of this mode are moving in the same direction and no nodes can be seen. 

The third and fourth mode on the right hand side of Fig. 4 appear similar to a plane wave oscillation of electrons, 
but trapped inside the cluster. To identify a wave number of the plane wave oscillation, a Fourier decomposition in 
z-dircction was done. A maximum at k = 1.6 nm _1 and k — 4.7 nm" 1 , respectively, appears, which belong to the 
wavelengths of the plane wave oscillations. The second resonance structure of Fig. 4, looks like a mix of the first and 
the third resonance structure. The plane wave oscillation with higher wavenumber was only found in the large cluster 
with 1000 ions. All other modes can be seen in smaller clusters as well. 

We want to point out one further feature of the mode spectra in Fig. 3 (right). The dashed red line represents 
two resonance structures with exactly the same eigenvalues at all frequencies. The eigenvectors are orthogonal. Both 
eigenvectors are characterized by the same spherical harmonic function Yi t i(JS, <p) but have a phase shift in (^-direction: 
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Y 1 \ '(6, (f) — Y^{9, 0+1). Further degenerations are obtained for lower eigenvalues. 

All eigenvectors ^(r, to) are decomposed into a superposition of spherical Bessel functions ji{k n ^r) with a set of 
ordinal numbers n. No leading ordinal number n was found, which characterizes the spatial resonance structure in r 
direction. 



A. Resonance frequency of the rigid oscillation 

The total current density ACF shown in Fig. 1 as well as the leading eigenvalue in Fig. 2 show a strong resonance 
at the frequency ss 4.42 fs _1 . This resonance belongs to the dipole-like mode with the eigenvector shown in Fig. 
4 on the left hand side. We will now analyze this collective excitation mode in terms of a rigid oscillation. 

The electrons with density profile n e {r) are assumed to move nearly rigidly in the external potential V ex t,ei(f) due 
to the fixed ions. The potential energy of the electrons due to a small shift with respect to the ions reads 



U e {z) = J d 3 fn c {f)V cxtM (f- 



ze z ). (22) 



The change of the potential energy U e (z) in z direction is due to the restoring force on the electron profile. In harmonic 
approximation of the equation of motion, the resonance frequency is identified as 

2 _ d 2 U c {z) 



m c N c uj^ — 



(23) 

z=0 



For small rigid shifts z — > 0, assuming radially dependent electron density profiles and external potentials in Eq.(22) 
the angular dependence of the potential energy calculation can be executed. The resonance frequency Eq.(23) is given 
according to 

"^dkr dr "= wr2 (^-» w+2 ^)^ < 24 » 

As a first example for a density profile, we assume a homogeneously charged ion sphere with radius R{ = 
(3iVi/(47rni)) 1 / 3 and an electron sphere with radius R c = (3-/V e /(47rn e )) 1/ ' 3 just on the density. The densities of 
the electron and ion spheres is equal (n c = rii). Therefore, the difference of ion and electron radius is determined by 
the cluster charge, basically the difference of the simulated electron number N e and ion number iV;. Thus, in the case 
of positive charged clusters, as discussed here, the electron sphere radius is smaller than the ion radius (R e < Ri). 
The error function potential Eq.(9) was taken as electron- ion-interaction potential for the calculation of the resonance 
frequency, as it was used for the MD simulation as well. The resonance frequency than reads 



Wp(i?i,i? c ) = u 



Mie 



Rf + i? 3 crf ( Rj + RA _ Rf-R 3 c ^ { (Ri - R c 



2i? 3 V A / 2iS V A 



y-A(i??+i? c 2 ) 



sinh I ^'f^ c ) — Ai?i_R c cosh 



A 2 / V A 2 



(25) 



In the limit of large clusters with high number of ions the resonance frequency equals the Mie frequency 
linifj.^oo cjft(i?j) = uJMic- In this case the sphere radii have nearly the same size (R e — > R\) and the system is 
nearly neutral. The limit for very small cluster, till down to just one atom, depends on the pseudopotential. In our 
case, the resonance frequency lim^^i wri(Ai) = if^ z^tx^m 1S ^ ue *° tne oscma ti° n 01 a single electron in the ionic 
error-function pseudo-potential Eq.(9). 

In Fig. 5 (left) the dependence of the resonance frequency wr of the dipole-like mode is shown in dependence on 
the size of the ion sphere. Empty circles represent results from MD simulations for a Nas5 cluster at n\ = 2.15 • 10 22 
cm~ 3 as well as for the Na 30 g cluster and Naiooo cluster at rii = 2.80- 10 22 cm~ 3 . The cluster size dependent resonance 
frequency was calculated using Eq.(25) for ion densities of n ; = 2.15 • 10 21 cm~ 3 (solid red line) and n\ — 2.80 • 10 21 
cm~ 3 (solid black line). The limit of both cases for large clusters, the Mie frequency a; 2 ^ = 3 e So ^ e , is given as well as 
dotted lines colored accordingly. 

Additionally, the electron density profile n e (r) was deducted from MD simulations for all cluster sizes and used to 
derive the resonance frequency wr solving Eq.(24) numerically As a result, one is able to identify the resonance fre- 
quency of the dipole-like mode with a deviation to the direct simulation results of less than 5%. Using homogeneously 
a charged ion and electron sphere leads to reasonable agreement in the limits of large clusters as well as for small 
clusters. However, the spatial structure of the density profile more accurate calculations of the resonance frequency 
in the intermediate cluster size regime. 
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FIG. 5. Left: Cluster size dependent resonance frequency ujB,(Ri)- Simulation results (empty circles) and analytical calculations 
(solid lines) using Eq.(25) are shown for ni = 2.15 • 10 22 cm -3 (red) and n\ = 2.80 • 10 22 cm~ 3 (black). The limit for large 
clusters is given as well (dashed line). Numerical calculations using Eq.(24) are presented (black and red dots). Right: Dispersion 
Eq.(32) of a plane wave in a homogeneously charged sphere (solid lines) for n e = 2.80 • 10 22 cm~ 3 as well as simulation results 
(empty symbols). 



B. Dispersion of the plane wave mode 



While the Mie-like resonance, discussed in the previous subsection is almost spherically symmetric, we obtain an 
increasing plane wave character of the dipole-like mode with increasing frequency. The third resonance frequency 
of the total current density ACF at wr = 8.14 fs -1 , see Fig. 1, is mainly caused by a plane wave like eigenvector 
as shown in Fig. 4 (right). Here, oscillations of electrons in opposite directions must be taken into account for the 
analytical calculation of the resonance frequency this. The electron motion is treated as a hydrodynamical liquid 
using the Euler equation 

d l(f,t) \->t-> *.\^->r-> ^ 1 a r->4\ n e(r,t) 



dt 



-div 



j(f, t) ® v(f, t) gradp(r, t) — grad V ext (f, t), (26) 



where j(r,t) — n c (r,t)v(r,t) is the spatially resolved current density of the electrons, p(r,t) is the pressure of the 
electron gas and V ext — K^xt.ei + ^ext.ee is the external potential, composed of contributions from the electrons and ions. 
First, one is able to linearize the Euler equation due to a small perturbation in z-direction restricted to longitudinal 
effects 

Using the following ansatz 

mt) = 8je z ^ k *-« t \ (27) 
v z (r,t) = 5ve z e^ kz - ut \ (28) 
n c (f,t) = n e , 0 (r) + *n 6 e i < fc *- wt >, (29) 
V ext (r, t) = Kxt,o + SV ext (r)e^ kz -^\ (30) 

for a small perturbation in ^-direction restricting ourselves to longitudinal effects, one is able to linearize the Euler 
equation. The system is assumed to be in LTE, described by the quantities n e ^(r), jo(r) =0, vo(r) =0 as well as 
V ext fi(r). One has to consider now electrons moving in an external field of ions and electrons to construct the external 
potential 

Kxt(r, t) = Kxt,o(r) + SV oxt {r, t). (31) 

Finally, the external potential has a equilibrium part and a perturbative part 6 V ext (r, t) , which is mainly dependent 
on the linear density perturbation 5n e (r, t). 

Assuming Boltzmann distribution to express the ideal gas pressure of the electrons p(r, t) via the electron density 
and using the equation of continuity, one is able to express the Euler equation in terms of linear perturbations of the 
density. Thus, the equilibrium part of the external potential compensates the pressure term on the right hand side 
of Eq.(26). Because we are interested in linear perturbations of the Euler equation, only the third term on the right 
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hand side of Eq.(26) remains, which is connected to the external potential. Finally, leading all terms of the Eulcr 
equation Eq.(26) back to a linear density fluctuation Sn e (r,t), one ends up with 

w 2(fc) = ^! e -!(fcA 2 +4iK e ) 

-e^ [l+c 2ifc ^]erf (^)) ■ (32) 

This relation leads to real solutions for the resonance frequencies only for standing waves with k n — nir/R e and scales 
with the plasma frequency co p i. 

From the eigenvector of the plane wave mode, one can derive the wave number k = ir/R e , which corresponds to 
n = 1. This means the dispersion of the plane wave mode is determined by the radius of the electron cloud. Results 
for this case are shown in Fig. 5 (right) for different cluster sizes and are compared with simulation data. For the 
cluster with 1000 ions a plane wave mode with k — 3ir/R e and n = 3, respectively, was found as well. The simulation 
data for the Naiooo cluster fit the dispersion Eq.(32) as well. Deviations of the plane wave resonance for smaller 
clusters are related to the radial dependence of the electron density profile. 

VI. CONCLUSION 

We have investigated collective excitation modes of a nano plasma in highly excited metal clusters. The collective 
excitation of electrons inside the cluster are obtained from bi-local current density correlation functions by solving the 
eigenvalue problem of the current density correlation matrix. Therefore, the local current density j(f,t) for excited 
clusters of 55 up to 1000 ions with densities of n; = 2.15 • 10 22 cm' 3 as well as 2.80 • 10 22 cm" 3 and temperatures 
of T e — 1 eV has been calculated using RMD simulations. Pseudo-potentials of sodium were used to calculate the 
electron dynamics without consideration of degeneration effects. This limits us to temperatures of a few T e > 1 eV. 
For the analysis of electron dynamics at temperatures less than T e = 1 eV, the inclusion of quantum effects for the 
calculation of the local current density j(f,t) of cluster electrons is an open question at this point. It is useful to go 
beyond classical description to discuss for example cold, non-excited clusters. 

The spectrum of a dipole-like mode was investigated. The leading mode more in detail. Using analytical calculations, 
it was possible to relate the position of resonance modes in the frequency domain to their spatial mode structure. 
The width of mode resonances and the role of collision-less damping effects as well as the collision frequency needs 
to be investigated in the future. The systematic change of the collision frequency with cluster size up to the bulk 
limit remains an interesting field. A first step was done here, presenting the cluster size dependence of resonance 
frequencies. The experimental verification of modes beyond the Mie resonance are desirable. Collective effects of 
electron motion have to be taken into account when analyzing ultraviolet (UPS) or x-ray photoelectron spectroscopy 
(XPS) experiments, as [22]. 
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